4.1.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- dat
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "bw_gain"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable     n  mean    sd
##   <chr>     <fct>    <dbl> <dbl> <dbl>
## 1 CTRL      bw_gain     12  17.1  2.26
## 2 PFOS      bw_gain     12  17.2  4.04
## 3 VAN       bw_gain     12  17.6  2.69
## 4 VAN+PFOS  bw_gain     12  17.5  2.65

4.1.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R01             1 no    no     310   321.  325.  339.  350    354
## 2 PFOS      R30            18 yes   no     258.  262.  265.  271.  270.   274
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.951  0.0457

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.530 0.664

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers, has equal variance and is normally distributed without the outliers according to Shapiro-Wilk test. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.1.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd     F     p p<.05   ges
## 1 treatment   3  44 0.064 0.979       0.004

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term   group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr>  <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>       
## 1 treat… CTRL   PFOS            0   0.111     -3.15      3.37 1     ns          
## 2 treat… CTRL   VAN             0   0.463     -2.79      3.72 0.981 ns          
## 3 treat… CTRL   VAN+P…          0   0.374     -2.88      3.63 0.99  ns          
## 4 treat… PFOS   VAN             0   0.353     -2.90      3.61 0.991 ns          
## 5 treat… PFOS   VAN+P…          0   0.264     -2.99      3.52 0.996 ns          
## 6 treat… VAN    VAN+P…          0  -0.0893    -3.35      3.17 1     ns

4.2.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- subset(dat, !is.na(cecum_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "cecum_norm"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 CTRL      cecum_norm    11 1     0.148
## 2 PFOS      cecum_norm    12 0.944 0.164
## 3 VAN       cecum_norm    12 1.88  0.24 
## 4 VAN+PFOS  cecum_norm    11 2.07  0.201

4.2.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 VAN       R15            27 no    yes    268.  277.  283.  290.  296.   300
## 2 VAN       R24            36 no    yes    281.  286.  294.  305.  309.   312
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.984   0.753

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    42     0.416 0.742

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised cecum weight data has two outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.2.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd       F        p p<.05   ges
## 1 treatment   3  42 106.226 1.21e-19     * 0.884

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1 group2   null.value estimate conf.low conf.high    p.adj
## * <chr>     <chr>  <chr>         <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
## 1 treatment CTRL   PFOS              0  -0.0562  -0.271      0.158 8.96e- 1
## 2 treatment CTRL   VAN               0   0.884    0.669      1.10  1.42e-12
## 3 treatment CTRL   VAN+PFOS          0   1.07     0.850      1.29  1.06e-12
## 4 treatment PFOS   VAN               0   0.940    0.730      1.15  1.09e-12
## 5 treatment PFOS   VAN+PFOS          0   1.13     0.911      1.34  1.06e-12
## 6 treatment VAN    VAN+PFOS          0   0.186   -0.0287     0.400 1.1 e- 1
## # ℹ 1 more variable: p.adj.signif <chr>

4.3.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove NA in the data column
dat.clean <- subset(dat, !is.na(liver_norm))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "liver_norm"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 CTRL      liver_norm    12 1     0.055
## 2 PFOS      liver_norm    12 1.08  0.055
## 3 VAN       liver_norm    12 0.934 0.044
## 4 VAN+PFOS  liver_norm    12 1.03  0.065

4.3.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Data contains zero outliers.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.985   0.778

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.430 0.733

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that normalised liver weight data has two outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

4.3.1.3 ANOVA One-Way test

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd      F        p p<.05  ges
## 1 treatment   3  44 15.909 3.79e-07     * 0.52

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1 group2   null.value estimate conf.low conf.high       p.adj
## * <chr>     <chr>  <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
## 1 treatment CTRL   PFOS              0   0.0846   0.0246   0.145   0.0027     
## 2 treatment CTRL   VAN               0  -0.0663  -0.126   -0.00623 0.0254     
## 3 treatment CTRL   VAN+PFOS          0   0.0352  -0.0249   0.0953  0.409      
## 4 treatment PFOS   VAN               0  -0.151   -0.211   -0.0909  0.000000182
## 5 treatment PFOS   VAN+PFOS          0  -0.0494  -0.110    0.0107  0.14       
## 6 treatment VAN    VAN+PFOS          0   0.102    0.0415   0.162   0.000269   
## # ℹ 1 more variable: p.adj.signif <chr>

4.3.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% difference",limits = c(0.75,1.5),breaks = seq(0.75,1.5,0.25), labels = function(x) paste0(x*100, "%")) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(1.35,1.4,1.45,1.5))
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Normalized liver weights at day 8

5.1.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove rows with NA
dat.clean <- subset(dat, !is.na(pfos_serum4_ugml))
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "pfos_serum4_ugml"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable             n   mean    sd
##   <chr>     <fct>            <dbl>  <dbl> <dbl>
## 1 CTRL      pfos_serum4_ugml    11  0     0.001
## 2 PFOS      pfos_serum4_ugml    12  9.17  2.01 
## 3 VAN       pfos_serum4_ugml    11  0.001 0.001
## 4 VAN+PFOS  pfos_serum4_ugml    10 10.0   1.18

5.1.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 5 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R06             6 no    no     274.  276.  285.  296.  297.   302
## 2 CTRL      R11            11 no    no     271.  278.  287.  294.  293.   299
## 3 PFOS      R29            17 yes   no     239.  243.  248.  255.  259.   267
## 4 VAN       R13            25 no    yes    218.  222.  228.  234   231.   241
## 5 VAN       R21            33 no    yes    262.  268.  274.  285.  281.   291
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic    p.value
##   <chr>                <dbl>      <dbl>
## 1 residuals(model)     0.807 0.00000412

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic        p
##   <int> <int>     <dbl>    <dbl>
## 1     3    40      7.48 0.000436

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers and has equal variance, however falls short on the Shapiro-Wilk test of normality and is therefore not normally distributed. Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.1.1.3 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.                  n statistic    df           p method        
## * <chr>            <int>     <dbl> <int>       <dbl> <chr>         
## 1 pfos_serum4_ugml    44      35.4     3 0.000000101 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.                  n effsize method  magnitude
## * <chr>            <int>   <dbl> <chr>   <ord>    
## 1 pfos_serum4_ugml    44   0.809 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
## * <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 pfos_serum4_… CTRL   PFOS      11    12    3.86   1.15e-4 1.85e-4 ***         
## 2 pfos_serum4_… CTRL   VAN       11    11    0.0172 9.86e-1 9.86e-1 ns          
## 3 pfos_serum4_… CTRL   VAN+P…    11    10    4.53   5.87e-6 1.91e-5 ****        
## 4 pfos_serum4_… PFOS   VAN       12    11   -3.84   1.23e-4 1.85e-4 ***         
## 5 pfos_serum4_… PFOS   VAN+P…    12    10    0.863  3.88e-1 4.66e-1 ns          
## 6 pfos_serum4_… VAN    VAN+P…    11    10    4.51   6.36e-6 1.91e-5 ****

5.1.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "ug/mL",limits = c(0,20),breaks = seq(0,20,5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(14,17,15,14))
p

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 16 rows containing missing values (`geom_point()`).

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
## Removed 16 rows containing missing values (`geom_point()`).

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 4 rows containing non-finite values (`stat_boxplot()`).
## Removed 16 rows containing missing values (`geom_point()`).

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS conc. in ug/mL at day 4

5.1.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum4_mg"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum4_mg))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.1.2.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable       statistic     p
##   <chr>     <chr>              <dbl> <dbl>
## 1 PFOS      pfos_serum4_mg     0.955 0.717
## 2 VAN+PFOS  pfos_serum4_mg     0.930 0.446

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20    0.0274 0.870

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.1.2.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1  -0.0171     0.165     0.183 pfos_s… PFOS   VAN+P…    12    10     -1.23 0.232
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

The output provides:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• statistic: Test statistic used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the estimated mean or difference in means depending on whether it was a one-sample test or a two-sample test.

• estimate1, estimate2: show the mean values of the two groups, respectively, for independent samples t-tests.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval.

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.            group1 group2   effsize    n1    n2 magnitude
## * <chr>          <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum4_mg PFOS   VAN+PFOS  -0.528    12    10 moderate

5.1.2.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "mg PFOS",limits = c(0,0.30),breaks = seq(0,0.30,0.1)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(0.28))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Total mg PFOS in blood volume

5.1.3.1 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 1 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 PFOS      R29            17 yes   no     239.  243.  248.  255.  259.   267
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable        statistic     p
##   <chr>     <chr>               <dbl> <dbl>
## 1 PFOS      pfos_serum4_pct     0.937 0.456
## 2 VAN+PFOS  pfos_serum4_pct     0.939 0.547

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20      1.43 0.246

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

5.1.3.2 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1   -0.626      6.75      7.37 pfos_s… PFOS   VAN+P…    12    10     -1.17 0.254
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

The output provides:

• .y.: the y variable used in the test.

• group1,group2: the compared groups in the pairwise tests.

• statistic: Test statistic used to compute the p-value.

• df: degrees of freedom.

• p: p-value.

• p.adj: the adjusted p-value.

• method: the statistical test used to compare groups.

• p.signif, p.adj.signif: the significance level of p-values and adjusted p-values, respectively.

• estimate: estimate of the effect size. It corresponds to the estimated mean or difference in means depending on whether it was a one-sample test or a two-sample test.

• estimate1, estimate2: show the mean values of the two groups, respectively, for independent samples t-tests.

• alternative: a character string describing the alternative hypothesis.

• conf.low,conf.high: Lower and upper bound on a confidence interval.

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.             group1 group2   effsize    n1    n2 magnitude
## * <chr>           <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum4_pct PFOS   VAN+PFOS  -0.503    12    10 moderate

5.1.3.3 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% of total dosed PFOS", limits = c(3,10),breaks = seq(3,10,1)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE) #, y.position = c(1.35,1.4,1.45,1.5))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum day 4 in pct. of total

5.2.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove rows with NA
dat.clean <- subset(dat, !is.na(pfos_serum8_ugml))
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "pfos_serum8_ugml"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable             n   mean     sd
##   <chr>     <fct>            <dbl>  <dbl>  <dbl>
## 1 CTRL      pfos_serum8_ugml    12  0.016  0.03 
## 2 PFOS      pfos_serum8_ugml    12 36.3   15.5  
## 3 VAN       pfos_serum8_ugml    12  0.011  0.021
## 4 VAN+PFOS  pfos_serum8_ugml    12 32.2   10.7

5.2.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 7 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R05             5 no    no     214   219.  222.  229.  231.   237
## 2 CTRL      R09             9 no    no     273.  278.  290.  290.  296.   302
## 3 CTRL      R12            12 no    no     297.  304.  316.  322.  324.   334
## 4 VAN       R14            26 no    yes    246.  256.  260.  267.  270.   274
## 5 VAN       R16            28 no    yes    256.  260   268.  275.  273    279
## 6 VAN       R24            36 no    yes    281.  286.  294.  305.  309.   312
## 7 VAN+PFOS  R47            47 yes   yes    242.  249.  255.  263.  267.   271
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic    p.value
##   <chr>                <dbl>      <dbl>
## 1 residuals(model)     0.823 0.00000461

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic         p
##   <int> <int>     <dbl>     <dbl>
## 1     3    44      8.92 0.0000993

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers and has equal variance, however falls short on the Shapiro-Wilk test of normality and is therefore not normally distributed. Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.2.1.3 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.                  n statistic    df            p method        
## * <chr>            <int>     <dbl> <int>        <dbl> <chr>         
## 1 pfos_serum8_ugml    48      37.3     3 0.0000000402 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.                  n effsize method  magnitude
## * <chr>            <int>   <dbl> <chr>   <ord>    
## 1 pfos_serum8_ugml    48   0.779 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.           group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
## * <chr>         <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 pfos_serum8_… CTRL   PFOS      12    12    4.37   1.26e-5 3.78e-5 ****        
## 2 pfos_serum8_… CTRL   VAN       12    12   -0.0899 9.28e-1 9.28e-1 ns          
## 3 pfos_serum8_… CTRL   VAN+P…    12    12    4.17   3.02e-5 4.53e-5 ****        
## 4 pfos_serum8_… PFOS   VAN       12    12   -4.46   8.32e-6 3.78e-5 ****        
## 5 pfos_serum8_… PFOS   VAN+P…    12    12   -0.195  8.46e-1 9.28e-1 ns          
## 6 pfos_serum8_… VAN    VAN+P…    12    12    4.26   2.03e-5 4.05e-5 ****

5.2.1.4 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "ug/mL",limits = c(0,80),breaks = seq(0,80,10)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = TRUE, y.position = c(72,80,75,72))
p

## Warning: Removed 11 rows containing missing values (`geom_point()`).

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 11 rows containing missing values (`geom_point()`).

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: Removed 11 rows containing missing values (`geom_point()`).

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS conc. in ug/mL at day 8

5.2.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum8_mg"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum8_mg))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.2.2.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable       statistic     p
##   <chr>     <chr>              <dbl> <dbl>
## 1 PFOS      pfos_serum8_mg     0.890 0.119
## 2 VAN+PFOS  pfos_serum8_mg     0.902 0.168

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    22      1.98 0.173

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.2.2.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1   0.0931     0.717     0.624 pfos_s… PFOS   VAN+P…    12    12     0.853 0.403
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.            group1 group2   effsize    n1    n2 magnitude
## * <chr>          <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum8_mg PFOS   VAN+PFOS   0.348    12    12 small

5.2.2.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "mg PFOS",limits = c(0,2),breaks = seq(0,2,0.5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(1.75))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

Total mg PFOS in blood volume

5.2.3.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_serum8_pct"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes" & !rat_name == "R47")

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_serum8_pct))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.2.3.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

##  [1] treatment         rat_name          ordering          pfos             
##  [5] van               bw_0              bw_1              bw_2             
##  [9] bw_3              bw_4              bw_5              bw_6             
## [13] bw_7              bw_8              bw_gain           cecum_wt         
## [17] cecum_wt_bw       cecum_norm        liver_wt          liver_wt_bw      
## [21] liver_norm        tot_pfos4         blood_vol4_mL     pfos_serum4_ugml 
## [25] pfos_serum4_ug    pfos_serum4_mg    pfos_serum4_pct   tot_pfos8        
## [29] blood_vol8_mL     pfos_serum8_ugml  pfos_serum8_ug    pfos_serum8_mg   
## [33] pfos_serum8_pct   pfos_change48_pct pfos_liver_ugg    pfos_liver_mg    
## [37] pfos_liver_pct    acetic            formic            propanoic        
## [41] m2_propanoic      butanoic          m3_butanoic       pentanoic        
## [45] m4_pentanoic      hexanoic          heptanoic         is.outlier       
## [49] is.extreme       
## <0 rækker> (eller 0-længde row.names)

Any extreme outliers can be bad samples or errors in data entry. If outliers compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable        statistic      p
##   <chr>     <chr>               <dbl>  <dbl>
## 1 PFOS      pfos_serum8_pct     0.867 0.0594
## 2 VAN+PFOS  pfos_serum8_pct     0.899 0.181

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    21      2.86 0.106

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

If the p-value of the Levene’s test is significant, it suggests that there is a significant difference between the variances of the two groups. In such case we should use Welch t-test, which doesn’t assume the equality of the two variances (var.equal=FALSE). If the Levene’s test is non-significant we can perform a Student t-test (var.equal=TRUE).

No outliers were identified. Data is normally distributed and has equal variance. Hence we use t-test.

5.2.3.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1     2.05      11.8      9.74 pfos_s… PFOS   VAN+P…    12    11      1.29 0.212
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.             group1 group2   effsize    n1    n2 magnitude
## * <chr>           <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_serum8_pct PFOS   VAN+PFOS   0.537    12    11 moderate

5.2.3.4 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% of total dosed PFOS", limits = c(5,25),breaks = seq(5,25,5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE, y.position = c(24))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 60, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum day 8 in pct. of total

5.3.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Set names of variables
PREDICTOR <- "treatment"
OUTCOME <- "pfos_change48_pct"
SUBJECT <- "rat_name"

# Subset to a specific varible
dat.clean <- subset(dat, pfos == "yes") # add following to subset() to remove the outliers: & !rat_name %in% c("R47","R27"))

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(pfos_change48_pct))

# Will yoou run a paired test? (set variable to `TRUE` or `FALSE`)
PAIRED <- FALSE

# Create formula
FORMULA <- as.formula(paste(OUTCOME, PREDICTOR, sep = "~"))

# Sort data for paired test
if (PAIRED) {
  # Order data
  dat.clean <- arrange(dat.clean, !!sym(SUBJECT))
  
  # Remove unpaired samples
  dat.clean <- dat.clean %>% 
    group_by(!!sym(SUBJECT)) %>%
    filter(n() != 1) %>%
    arrange(!!sym(PREDICTOR), !!sym(SUBJECT)) %>%
    droplevels() %>% 
    ungroup()
}

5.3.1.2 Assumptions and preliminary tests

The two-samples t-tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# identify outliers
dat.clean %>%
  group_by(!!sym(PREDICTOR)) %>%
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 PFOS      R27            15 yes   no     270.  280.  284.  291.  290.   296
## 2 VAN+PFOS  R47            47 yes   yes    242.  249.  255.  263.  267.   271
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Any extreme outliers can be bad samples or errors in data entry. If outliers, compare a test with and without the outlier to determine if it is important, or run a non-parametric Wilcoxon test.

Check normality by groups
The normality assumption can be checked by computing the Shapiro-Wilk test for each group. If the data is normally distributed, the p-value should be greater than 0.05. You can also create QQ plots for each group. QQ plot draws the correlation between a given data and the normal distribution.

If your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.

Consequently, we should not rely on only one approach for assessing the normality. A better strategy is to combine visual inspection and statistical test.

# Run Shapiro test
dat.clean %>% 
  group_by(!!sym(PREDICTOR)) %>%
  shapiro_test(!!sym(OUTCOME))

## # A tibble: 2 × 4
##   treatment variable          statistic      p
##   <chr>     <chr>                 <dbl>  <dbl>
## 1 PFOS      pfos_change48_pct     0.894 0.132 
## 2 VAN+PFOS  pfos_change48_pct     0.805 0.0167

# Create QQplot
ggqqplot(dat.clean, x = OUTCOME, facet.by = PREDICTOR)

If both Shapiro test has p > 0.05 and/ or the QQplot follows the reference line the data follows a normal distribution.

If the data does not follow the normal distribution run a Wilcoxon Rank-sum test

Check the equality of variances
This can be done using the Levene’s test. If the variances of groups are equal, the p-value should be greater than 0.05.

# Run test
dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     1    20     0.639 0.434

# Save output
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

Two outliers were identified (sample for R27 and R47). Analysis result and test method is similar with and without outliers. Data is normally distributed and has equal variance. Hence we use t-test.

5.3.1.3 PERFORM TEST

T-test
We are now ready to perform the test

stat.test <- dat.clean %>% 
  t_test(FORMULA,
         var.equal = EQUAL.VAR,
         detailed = TRUE,
         paired = FALSE,
         alternative = "two.sided") %>%
  add_significance()
stat.test

## # A tibble: 1 × 16
##   estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic     p
##      <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl>
## 1     85.1      341.      256. pfos_c… PFOS   VAN+P…    12    10      1.11 0.281
## # ℹ 6 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>

Effect size
The effect size is calculated as Cohen’s D

dat.clean %>% cohens_d(FORMULA, 
                       var.equal = EQUAL.VAR,
                       paired = FALSE)

## # A tibble: 1 × 7
##   .y.               group1 group2   effsize    n1    n2 magnitude
## * <chr>             <chr>  <chr>      <dbl> <int> <int> <ord>    
## 1 pfos_change48_pct PFOS   VAN+PFOS   0.475    12    10 small

5.3.1.4 Conclusion

5.3.1.5 Create figure

# Prepare stats
stat.test <- stat.test %>% add_xy_position(x = PREDICTOR)

# Create point plot with mean and SD
data_summary <- function(x) {
  m <- mean(x)
  ymin <- m-sd(x)
  ymax <- m+sd(x)
  return(c(y=m,ymin=ymin,ymax=ymax))
}
data_summary_collapsed <- function(x) {
  m <- mean(x)
  ymin <- m
  ymax <- m
  return(c(y=m,ymin=ymin,ymax=ymax))
}

p <- ggplot(dat.clean, aes(x = .data[[PREDICTOR]], y = .data[[OUTCOME]], color = .data[[PREDICTOR]])) +
  stat_summary(fun.data = data_summary_collapsed, geom = "crossbar", color = "black", width = 0.5, linewidth = 0.3) +
  stat_summary(fun.data = data_summary, geom = "errorbar", color = "black", width = 0.15, linewidth = 0.5) +
  geom_point(position = position_jitterdodge(dodge.width = 0.6, jitter.width = 0.4), size = 2, colour = "black", shape = 21, stroke = 0.5, aes(fill = treatment)) +
  scale_fill_manual(values = params$COL) +
  scale_y_continuous(name = "% change", limits = c(100,900),breaks = seq(100,900,100), labels = function(x) paste0(x, "%")) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment") +
  theme_pubr()
p

# Alternative: Create boxplot
# p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
#           fill = PREDICTOR,
#           add =  "jitter",
#           add.params = list(size = 1)) +
#   scale_fill_manual(values = params$COL) +
#   scale_y_continuous(name = "% change", limits = c(100,900),breaks = seq(100,900,100)) +
#   labs(fill = "Treatment") +
#   scale_x_discrete(name = "Treatment")

p <- p + stat_pvalue_manual(stat.test, tip.length = 0, hide.ns = FALSE) #, y.position = c(1.35,1.4,1.45,1.5))
p2 <- p + labs(subtitle = get_test_label(stat.test, detailed = TRUE))
p2

p

# Plot for saving without legend
p3 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 90, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/",OUTCOME,"_plot.pdf"), p3, device = "pdf", dpi = 300, units = "mm", width = 70, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS serum level change from day 4 to 8 in pct.

5.3.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Color scheme
COL <- c("#61d46b","#ffe900","#31b44b","#efc000")

# Subset data
dat.sub <- subset(dat, pfos == "yes")

# Create data frame for data representation
dat.clean <- dat.sub %>% select(rat_name, treatment, pfos_serum4_ugml, pfos_serum8_ugml) %>%
  pivot_longer(., cols = c(pfos_serum4_ugml, pfos_serum8_ugml), names_to = "data_group", values_to = "conc")

# Create column for day of sampling
dat.clean <- transform(dat.clean, "day" = ifelse(dat.clean$data_group == "pfos_serum8_ugml","d8","d4"))

# Create ID column for easier handling
for (i in dat.sub$rat_name) {
  dat.clean$ID <- paste(dat.clean$day,"_",dat.clean$treatment)
}

# Order dataframe for analysis
dat.clean <- dat.clean[order(dat.clean$day),]

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(conc))
dat.clean

##    rat_name treatment       data_group  conc day            ID
## 1       R25      PFOS pfos_serum4_ugml  8.56  d4     d4 _ PFOS
## 3       R26      PFOS pfos_serum4_ugml  9.60  d4     d4 _ PFOS
## 5       R27      PFOS pfos_serum4_ugml  7.92  d4     d4 _ PFOS
## 7       R28      PFOS pfos_serum4_ugml  8.64  d4     d4 _ PFOS
## 9       R29      PFOS pfos_serum4_ugml 12.96  d4     d4 _ PFOS
## 11      R30      PFOS pfos_serum4_ugml  8.68  d4     d4 _ PFOS
## 13      R31      PFOS pfos_serum4_ugml  5.72  d4     d4 _ PFOS
## 15      R32      PFOS pfos_serum4_ugml  7.56  d4     d4 _ PFOS
## 17      R33      PFOS pfos_serum4_ugml  8.24  d4     d4 _ PFOS
## 19      R34      PFOS pfos_serum4_ugml 11.36  d4     d4 _ PFOS
## 21      R35      PFOS pfos_serum4_ugml  9.00  d4     d4 _ PFOS
## 23      R36      PFOS pfos_serum4_ugml 11.84  d4     d4 _ PFOS
## 25      R37  VAN+PFOS pfos_serum4_ugml  9.36  d4 d4 _ VAN+PFOS
## 27      R38  VAN+PFOS pfos_serum4_ugml 11.76  d4 d4 _ VAN+PFOS
## 29      R39  VAN+PFOS pfos_serum4_ugml 10.64  d4 d4 _ VAN+PFOS
## 31      R40  VAN+PFOS pfos_serum4_ugml 12.12  d4 d4 _ VAN+PFOS
## 33      R41  VAN+PFOS pfos_serum4_ugml  9.12  d4 d4 _ VAN+PFOS
## 35      R42  VAN+PFOS pfos_serum4_ugml  9.80  d4 d4 _ VAN+PFOS
## 37      R43  VAN+PFOS pfos_serum4_ugml 10.08  d4 d4 _ VAN+PFOS
## 43      R46  VAN+PFOS pfos_serum4_ugml  9.88  d4 d4 _ VAN+PFOS
## 45      R47  VAN+PFOS pfos_serum4_ugml  8.80  d4 d4 _ VAN+PFOS
## 47      R48  VAN+PFOS pfos_serum4_ugml  8.60  d4 d4 _ VAN+PFOS
## 2       R25      PFOS pfos_serum8_ugml 43.92  d8     d8 _ PFOS
## 4       R26      PFOS pfos_serum8_ugml 45.60  d8     d8 _ PFOS
## 6       R27      PFOS pfos_serum8_ugml 69.92  d8     d8 _ PFOS
## 8       R28      PFOS pfos_serum8_ugml 21.92  d8     d8 _ PFOS
## 10      R29      PFOS pfos_serum8_ugml 26.56  d8     d8 _ PFOS
## 12      R30      PFOS pfos_serum8_ugml 47.84  d8     d8 _ PFOS
## 14      R31      PFOS pfos_serum8_ugml 21.84  d8     d8 _ PFOS
## 16      R32      PFOS pfos_serum8_ugml 29.36  d8     d8 _ PFOS
## 18      R33      PFOS pfos_serum8_ugml 22.08  d8     d8 _ PFOS
## 20      R34      PFOS pfos_serum8_ugml 52.48  d8     d8 _ PFOS
## 22      R35      PFOS pfos_serum8_ugml 29.84  d8     d8 _ PFOS
## 24      R36      PFOS pfos_serum8_ugml 23.76  d8     d8 _ PFOS
## 26      R37  VAN+PFOS pfos_serum8_ugml 36.08  d8 d8 _ VAN+PFOS
## 28      R38  VAN+PFOS pfos_serum8_ugml 31.28  d8 d8 _ VAN+PFOS
## 30      R39  VAN+PFOS pfos_serum8_ugml 25.92  d8 d8 _ VAN+PFOS
## 32      R40  VAN+PFOS pfos_serum8_ugml 23.92  d8 d8 _ VAN+PFOS
## 34      R41  VAN+PFOS pfos_serum8_ugml 21.68  d8 d8 _ VAN+PFOS
## 36      R42  VAN+PFOS pfos_serum8_ugml 40.96  d8 d8 _ VAN+PFOS
## 38      R43  VAN+PFOS pfos_serum8_ugml 25.60  d8 d8 _ VAN+PFOS
## 40      R44  VAN+PFOS pfos_serum8_ugml 37.44  d8 d8 _ VAN+PFOS
## 42      R45  VAN+PFOS pfos_serum8_ugml 25.36  d8 d8 _ VAN+PFOS
## 44      R46  VAN+PFOS pfos_serum8_ugml 34.72  d8 d8 _ VAN+PFOS
## 46      R47  VAN+PFOS pfos_serum8_ugml 59.52  d8 d8 _ VAN+PFOS
## 48      R48  VAN+PFOS pfos_serum8_ugml 23.76  d8 d8 _ VAN+PFOS

# Set names of variables
PREDICTOR <- "ID"
OUTCOME <- "conc"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   ID            variable     n  mean    sd
##   <chr>         <fct>    <dbl> <dbl> <dbl>
## 1 d4 _ PFOS     conc        12  9.17  2.01
## 2 d4 _ VAN+PFOS conc        10 10.0   1.18
## 3 d8 _ PFOS     conc        12 36.3  15.5 
## 4 d8 _ VAN+PFOS conc        12 32.2  10.7

5.3.2.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = COL)
bxp

5.3.2.3 Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 8
##   ID            rat_name treatment data_group   conc day   is.outlier is.extreme
##   <chr>         <chr>    <chr>     <chr>       <dbl> <chr> <lgl>      <lgl>     
## 1 d4 _ PFOS     R29      PFOS      pfos_serum…  13.0 d4    TRUE       FALSE     
## 2 d8 _ VAN+PFOS R47      VAN+PFOS  pfos_serum…  59.5 d8    TRUE       FALSE

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.875 0.000151

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic       p
##   <int> <int>     <dbl>   <dbl>
## 1     3    42      6.46 0.00108

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that PFOS concentrations from day 4 and 8 has two outliers, has unequal variance, and falls short on the Shapiro-Wilk test of normality (not normally distributed). Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.3.2.4 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.       n statistic    df           p method        
## * <chr> <int>     <dbl> <int>       <dbl> <chr>         
## 1 conc     46      34.4     3 0.000000165 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.       n effsize method  magnitude
## * <chr> <int>   <dbl> <chr>   <ord>    
## 1 conc     46   0.747 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.   group1        group2     n1    n2 statistic       p   p.adj p.adj.signif
## * <chr> <chr>         <chr>   <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 conc  d4 _ PFOS     d4 _ V…    12    10     0.798 4.25e-1 5.10e-1 ns          
## 2 conc  d4 _ PFOS     d8 _ P…    12    12     4.68  2.92e-6 1.75e-5 ****        
## 3 conc  d4 _ PFOS     d8 _ V…    12    12     4.48  7.51e-6 2.25e-5 ****        
## 4 conc  d4 _ VAN+PFOS d8 _ P…    10    12     3.66  2.51e-4 5.02e-4 ***         
## 5 conc  d4 _ VAN+PFOS d8 _ V…    10    12     3.47  5.15e-4 7.73e-4 ***         
## 6 conc  d8 _ PFOS     d8 _ V…    12    12    -0.198 8.43e-1 8.43e-1 ns

5.3.2.5 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = COL,labels = c("PFOS day 4","VAN+PFOS day 4","PFOS day 8","VAN+PFOS day 8")) +
  scale_y_continuous(name = "ug/mL",limits = c(0,85),breaks = seq(0,85,10)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment", labels = c("PFOS\nday 4","VAN+PFOS\nday 4","PFOS\nday 8","VAN+PFOS\nday 8"))

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = FALSE, y.position = c(75,85,70,80))
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/pfos_day48_",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/pfos_day48_",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS conc. in ug/mL at day 8

5.3.3.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Color scheme
COL <- c("#61d46b","#ffe900","#31b44b","#efc000")

# Subset data
dat.sub <- subset(dat, pfos == "yes")

# Create data frame for data representation
dat.clean <- dat.sub %>% select(rat_name, treatment, pfos_serum4_mg, pfos_serum8_mg) %>%
  pivot_longer(., cols = c(pfos_serum4_mg, pfos_serum8_mg), names_to = "data_group", values_to = "mg")

# Create column for day of sampling
dat.clean <- transform(dat.clean, "day" = ifelse(dat.clean$data_group == "pfos_serum8_mg","d8","d4"))

# Create ID column for easier handling
for (i in dat.sub$rat_name) {
  dat.clean$ID <- paste(dat.clean$day,"_",dat.clean$treatment)
}

# Order dataframe for analysis
dat.clean <- dat.clean[order(dat.clean$day),]

# Remove rows with NA
dat.clean <- subset(dat.clean, !is.na(mg))
dat.clean

##    rat_name treatment     data_group        mg day            ID
## 1       R25      PFOS pfos_serum4_mg 0.1904292  d4     d4 _ PFOS
## 3       R26      PFOS pfos_serum4_mg 0.1592525  d4     d4 _ PFOS
## 5       R27      PFOS pfos_serum4_mg 0.1472486  d4     d4 _ PFOS
## 7       R28      PFOS pfos_serum4_mg 0.1480827  d4     d4 _ PFOS
## 9       R29      PFOS pfos_serum4_mg 0.2149079  d4     d4 _ PFOS
## 11      R30      PFOS pfos_serum4_mg 0.1498237  d4     d4 _ PFOS
## 13      R31      PFOS pfos_serum4_mg 0.1123500  d4     d4 _ PFOS
## 15      R32      PFOS pfos_serum4_mg 0.1359590  d4     d4 _ PFOS
## 17      R33      PFOS pfos_serum4_mg 0.1546747  d4     d4 _ PFOS
## 19      R34      PFOS pfos_serum4_mg 0.2121503  d4     d4 _ PFOS
## 21      R35      PFOS pfos_serum4_mg 0.1792512  d4     d4 _ PFOS
## 23      R36      PFOS pfos_serum4_mg 0.1811804  d4     d4 _ PFOS
## 25      R37  VAN+PFOS pfos_serum4_mg 0.1927112  d4 d4 _ VAN+PFOS
## 27      R38  VAN+PFOS pfos_serum4_mg 0.2344474  d4 d4 _ VAN+PFOS
## 29      R39  VAN+PFOS pfos_serum4_mg 0.1721467  d4 d4 _ VAN+PFOS
## 31      R40  VAN+PFOS pfos_serum4_mg 0.2381338  d4 d4 _ VAN+PFOS
## 33      R41  VAN+PFOS pfos_serum4_mg 0.1657651  d4 d4 _ VAN+PFOS
## 35      R42  VAN+PFOS pfos_serum4_mg 0.1652672  d4 d4 _ VAN+PFOS
## 37      R43  VAN+PFOS pfos_serum4_mg 0.2037289  d4 d4 _ VAN+PFOS
## 43      R46  VAN+PFOS pfos_serum4_mg 0.1702838  d4 d4 _ VAN+PFOS
## 45      R47  VAN+PFOS pfos_serum4_mg 0.1501491  d4 d4 _ VAN+PFOS
## 47      R48  VAN+PFOS pfos_serum4_mg 0.1332518  d4 d4 _ VAN+PFOS
## 2       R25      PFOS pfos_serum8_mg 1.0930000  d8     d8 _ PFOS
## 4       R26      PFOS pfos_serum8_mg 0.8230000  d8     d8 _ PFOS
## 6       R27      PFOS pfos_serum8_mg 1.3690000  d8     d8 _ PFOS
## 8       R28      PFOS pfos_serum8_mg 0.3980000  d8     d8 _ PFOS
## 10      R29      PFOS pfos_serum8_mg 0.4810000  d8     d8 _ PFOS
## 12      R30      PFOS pfos_serum8_mg 0.8450000  d8     d8 _ PFOS
## 14      R31      PFOS pfos_serum8_mg 0.4720000  d8     d8 _ PFOS
## 16      R32      PFOS pfos_serum8_mg 0.5750000  d8     d8 _ PFOS
## 18      R33      PFOS pfos_serum8_mg 0.4550000  d8     d8 _ PFOS
## 20      R34      PFOS pfos_serum8_mg 1.0550000  d8     d8 _ PFOS
## 22      R35      PFOS pfos_serum8_mg 0.6510000  d8     d8 _ PFOS
## 24      R36      PFOS pfos_serum8_mg 0.3890000  d8     d8 _ PFOS
## 26      R37  VAN+PFOS pfos_serum8_mg 0.7850000  d8 d8 _ VAN+PFOS
## 28      R38  VAN+PFOS pfos_serum8_mg 0.6670000  d8 d8 _ VAN+PFOS
## 30      R39  VAN+PFOS pfos_serum8_mg 0.4530000  d8 d8 _ VAN+PFOS
## 32      R40  VAN+PFOS pfos_serum8_mg 0.5010000  d8 d8 _ VAN+PFOS
## 34      R41  VAN+PFOS pfos_serum8_mg 0.4270000  d8 d8 _ VAN+PFOS
## 36      R42  VAN+PFOS pfos_serum8_mg 0.7600000  d8 d8 _ VAN+PFOS
## 38      R43  VAN+PFOS pfos_serum8_mg 0.5540000  d8 d8 _ VAN+PFOS
## 40      R44  VAN+PFOS pfos_serum8_mg 0.7520000  d8 d8 _ VAN+PFOS
## 42      R45  VAN+PFOS pfos_serum8_mg 0.4590000  d8 d8 _ VAN+PFOS
## 44      R46  VAN+PFOS pfos_serum8_mg 0.6440000  d8 d8 _ VAN+PFOS
## 46      R47  VAN+PFOS pfos_serum8_mg 1.0890000  d8 d8 _ VAN+PFOS
## 48      R48  VAN+PFOS pfos_serum8_mg 0.3980000  d8 d8 _ VAN+PFOS

# Set names of variables
PREDICTOR <- "ID"
OUTCOME <- "mg"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   ID            variable     n  mean    sd
##   <chr>         <fct>    <dbl> <dbl> <dbl>
## 1 d4 _ PFOS     mg          12 0.165 0.031
## 2 d4 _ VAN+PFOS mg          10 0.183 0.034
## 3 d8 _ PFOS     mg          12 0.717 0.32 
## 4 d8 _ VAN+PFOS mg          12 0.624 0.201

5.3.3.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = COL)
bxp

5.3.3.3 Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## [1] ID         rat_name   treatment  data_group mg         day        is.outlier
## [8] is.extreme
## <0 rækker> (eller 0-længde row.names)

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.896 0.000640

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic         p
##   <int> <int>     <dbl>     <dbl>
## 1     3    42      9.67 0.0000564

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that PFOS concentrations from day 4 and 8 has two outliers, has unequal variance, and falls short on the Shapiro-Wilk test of normality (not normally distributed). Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.3.3.4 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.       n statistic    df           p method        
## * <chr> <int>     <dbl> <int>       <dbl> <chr>         
## 1 mg       46      34.1     3 0.000000187 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.       n effsize method  magnitude
## * <chr> <int>   <dbl> <chr>   <ord>    
## 1 mg       46   0.741 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc

## # A tibble: 6 × 9
##   .y.   group1        group2     n1    n2 statistic       p   p.adj p.adj.signif
## * <chr> <chr>         <chr>   <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 mg    d4 _ PFOS     d4 _ V…    12    10     0.574 5.66e-1 6.79e-1 ns          
## 2 mg    d4 _ PFOS     d8 _ P…    12    12     4.62  3.92e-6 2.35e-5 ****        
## 3 mg    d4 _ PFOS     d8 _ V…    12    12     4.33  1.51e-5 4.54e-5 ****        
## 4 mg    d4 _ VAN+PFOS d8 _ P…    10    12     3.83  1.30e-4 2.60e-4 ***         
## 5 mg    d4 _ VAN+PFOS d8 _ V…    10    12     3.55  3.84e-4 5.75e-4 ***         
## 6 mg    d8 _ PFOS     d8 _ V…    12    12    -0.289 7.73e-1 7.73e-1 ns

5.3.3.5 Create figure

## Prepare statistical information:
pwc.adj <- pwc %>% 
  add_x_position(x = PREDICTOR) %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)

# Format for ggplot
if (sum(pwc.adj$p.adj.signif != "ns") == 0) {
  stat.sig <- pwc.adj %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
} else {
  stat.sig <- pwc.adj[pwc.adj$p.adj.signif != "ns",] %>%
    add_y_position(step.increase = 0.25) %>%
    mutate(y.position = seq(min(y.position), max(y.position),length.out = n()))
}

# Create plot
p <- ggboxplot(dat.clean, x = PREDICTOR, y = OUTCOME,
          fill = PREDICTOR,
          add =  "jitter",
          add.params = list(size = 1)) +
  scale_fill_manual(values = COL,labels = c("PFOS day 4","VAN+PFOS day 4","PFOS day 8","VAN+PFOS day 8")) +
  scale_y_continuous(name = "mg",limits = c(0,1.75),breaks = seq(0,1.75,0.5)) +
  labs(fill = "Treatment") +
  scale_x_discrete(name = "Treatment", labels = c("PFOS\nday 4","VAN+PFOS\nday 4","PFOS\nday 8","VAN+PFOS\nday 8")) +
  theme(axis.title.x = element_blank())

p <- p + stat_pvalue_manual(stat.sig, label = "p.adj.format",tip.length = 0, hide.ns = FALSE, y.position = c(1.635,1.75,1.4,1.52))
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

ggsave(filename = paste0("plots/animal_data/pfos/pfos_day48_",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/pfos/pfos_day48_",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

PFOS amount in mg at day 4 and 8

5.4.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
# Remove rows with NA
dat.clean <- subset(dat, !is.na(pfos_liver_ugg))
#dat.clean <- dat %>% select_if(~ !any(is.na(.)))
#dat.clean <- subset(dat, !dat$rat_name %in% c("R01","R30"))

# Set names of variables
PREDICTOR <- "treatment"#c("treatment","pfos","van")
OUTCOME <- "pfos_liver_ugg"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   treatment variable           n    mean     sd
##   <chr>     <fct>          <dbl>   <dbl>  <dbl>
## 1 CTRL      pfos_liver_ugg    12   0.199  0.173
## 2 PFOS      pfos_liver_ugg    12 176.    21.7  
## 3 VAN       pfos_liver_ugg    12   0.205  0.22 
## 4 VAN+PFOS  pfos_liver_ugg    12 196.    18.5

5.4.1.2 Visualise

Create a boxplot of the data.

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

#### Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 49
##   treatment rat_name ordering pfos  van    bw_0  bw_1  bw_2  bw_3  bw_4  bw_5
##   <chr>     <chr>       <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 CTRL      R11            11 no    no     271.  278.  287.  294.  293.   299
## 2 VAN       R16            28 no    yes    256.  260   268.  275.  273    279
## # ℹ 38 more variables: bw_6 <int>, bw_7 <dbl>, bw_8 <int>, bw_gain <dbl>,
## #   cecum_wt <dbl>, cecum_wt_bw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wt_bw <dbl>, liver_norm <dbl>, tot_pfos4 <dbl>, blood_vol4_mL <dbl>,
## #   pfos_serum4_ugml <dbl>, pfos_serum4_ug <dbl>, pfos_serum4_mg <dbl>,
## #   pfos_serum4_pct <dbl>, tot_pfos8 <dbl>, blood_vol8_mL <dbl>,
## #   pfos_serum8_ugml <dbl>, pfos_serum8_ug <dbl>, pfos_serum8_mg <dbl>,
## #   pfos_serum8_pct <dbl>, pfos_change48_pct <dbl>, pfos_liver_ugg <dbl>, …

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic  p.value
##   <chr>                <dbl>    <dbl>
## 1 residuals(model)     0.877 0.000119

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic          p
##   <int> <int>     <dbl>      <dbl>
## 1     3    44      12.2 0.00000631

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers and has equal variance, however falls short on the Shapiro-Wilk test of normality and is therefore not normally distributed. Therefore we use a non-parametric Kruskal-Wallis test with Dunn’s p-value adjustment.

5.4.1.3 Kruskal-Wallis test

res.aov <- dat.clean %>% kruskal_test(FORMULA)
res.aov

## # A tibble: 1 × 6
##   .y.                n statistic    df            p method        
## * <chr>          <int>     <dbl> <int>        <dbl> <chr>         
## 1 pfos_liver_ugg    48      36.4     3 0.0000000608 Kruskal-Wallis

dat.clean %>% kruskal_effsize(FORMULA)

## # A tibble: 1 × 5
##   .y.                n effsize method  magnitude
## * <chr>          <int>   <dbl> <chr>   <ord>    
## 1 pfos_liver_ugg    48   0.760 eta2[H] large

# pairwise comparisons
pwc <- dat.clean %>% 
  dunn_test(FORMULA, p.adjust.method = "fdr") 
pwc